home *** CD-ROM | disk | FTP | other *** search
/ MPEG Toolkit / MPEG Toolkit.iso / dos / mpegstat / jrevdct.c < prev    next >
Encoding:
C/C++ Source or Header  |  1997-01-01  |  46.5 KB  |  1,504 lines
  1. /* MPEGSTAT - analyzing tool for MPEG-I video streams
  2.  * 
  3.  * Technical University of Berlin, Germany, Dept. of Computer Science
  4.  * Tom Pfeifer - Multimedia systems project - pfeifer@fokus.gmd.de
  5.  *
  6.  * Jens Brettin, Harald Masche, Alexander Schulze, Dirk Schubert
  7.  *
  8.  * This program uses parts of the source code of the Berkeley MPEG player
  9.  *
  10.  * ---------------------------
  11.  *
  12.  * Copyright (c) 1993 Technical University of Berlin, Germany
  13.  *
  14.  * for the parts of the Berkeley player used:
  15.  *
  16.  * Copyright (c) 1992 The Regents of the University of California.
  17.  * All rights reserved.
  18.  *
  19.  * ---------------------------
  20.  *
  21.  * Permission to use, copy, modify, and distribute this software and its
  22.  * documentation for any purpose, without fee, and without written agreement is
  23.  * hereby granted, provided that the above copyright notices and the following
  24.  * two paragraphs appear in all copies of this software.
  25.  * 
  26.  * IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA 
  27.  * or the Technical University of Berlin BE LIABLE TO ANY PARTY FOR
  28.  * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT
  29.  * OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF
  30.  * CALIFORNIA or the Technical University of Berlin HAS BEEN ADVISED OF THE 
  31.  * POSSIBILITY OF SUCH DAMAGE.
  32.  * 
  33.  * THE UNIVERSITY OF CALIFORNIA and the Technical University of Berlin 
  34.  * SPECIFICALLY DISCLAIM ANY WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
  35.  * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
  36.  * PURPOSE.  THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE 
  37.  * UNIVERSITY OF CALIFORNIA and the Technical University of Berlin HAVE NO 
  38.  * OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, 
  39.  * OR MODIFICATIONS.
  40.  */
  41.  
  42. /*
  43.  * jrevdct.c
  44.  *
  45.  * Copyright (C) 1991, 1992, Thomas G. Lane.
  46.  * This file is part of the Independent JPEG Group's software.
  47.  * For conditions of distribution and use, see the accompanying README file.
  48.  *
  49.  * This file contains the basic inverse-DCT transformation subroutine.
  50.  *
  51.  * This implementation is based on an algorithm described in
  52.  *   C. Loeffler, A. Ligtenberg and G. Moschytz, "Practical Fast 1-D DCT
  53.  *   Algorithms with 11 Multiplications", Proc. Int'l. Conf. on Acoustics,
  54.  *   Speech, and Signal Processing 1989 (ICASSP '89), pp. 988-991.
  55.  * The primary algorithm described there uses 11 multiplies and 29 adds.
  56.  * We use their alternate method with 12 multiplies and 32 adds.
  57.  * The advantage of this method is that no data path contains more than one
  58.  * multiplication; this allows a very simple and accurate implementation in
  59.  * scaled fixed-point arithmetic, with a minimal number of shifts.
  60.  * 
  61.  * I've made lots of modifications to attempt to take advantage of the
  62.  * sparse nature of the DCT matrices we're getting.  Although the logic
  63.  * is cumbersome, it's straightforward and the resulting code is much
  64.  * faster.
  65.  *
  66.  * A better way to do this would be to pass in the DCT block as a sparse
  67.  * matrix, perhaps with the difference cases encoded.
  68.  */
  69.  
  70. #include <string.h>
  71. #include "video.h"
  72. #include "proto.h"
  73. #include "dither.h"
  74.  
  75. #define GLOBAL            /* a function referenced thru EXTERNs */
  76.   
  77. /* We assume that right shift corresponds to signed division by 2 with
  78.  * rounding towards minus infinity.  This is correct for typical "arithmetic
  79.  * shift" instructions that shift in copies of the sign bit.  But some
  80.  * C compilers implement >> with an unsigned shift.  For these machines you
  81.  * must define RIGHT_SHIFT_IS_UNSIGNED.
  82.  * RIGHT_SHIFT provides a proper signed right shift of an INT32 quantity.
  83.  * It is only applied with constant shift counts.  SHIFT_TEMPS must be
  84.  * included in the variables of any routine using RIGHT_SHIFT.
  85.  */
  86.   
  87. #ifdef RIGHT_SHIFT_IS_UNSIGNED
  88. #define SHIFT_TEMPS    INT32 shift_temp;
  89. #define RIGHT_SHIFT(x,shft)  \
  90.     ((shift_temp = (x)) < 0 ? \
  91.      (shift_temp >> (shft)) | ((~((INT32) 0)) << (32-(shft))) : \
  92.      (shift_temp >> (shft)))
  93. #else
  94. #define SHIFT_TEMPS
  95. #define RIGHT_SHIFT(x,shft)    ((x) >> (shft))
  96. #endif
  97.  
  98. /*
  99.  * This routine is specialized to the case DCTSIZE = 8.
  100.  */
  101.  
  102. #if DCTSIZE != 8
  103.   Sorry, this code only copes with 8x8 DCTs. /* deliberate syntax err */
  104. #endif
  105.  
  106.  
  107. /*
  108.  * A 2-D IDCT can be done by 1-D IDCT on each row followed by 1-D IDCT
  109.  * on each column.  Direct algorithms are also available, but they are
  110.  * much more complex and seem not to be any faster when reduced to code.
  111.  *
  112.  * The poop on this scaling stuff is as follows:
  113.  *
  114.  * Each 1-D IDCT step produces outputs which are a factor of sqrt(N)
  115.  * larger than the true IDCT outputs.  The final outputs are therefore
  116.  * a factor of N larger than desired; since N=8 this can be cured by
  117.  * a simple right shift at the end of the algorithm.  The advantage of
  118.  * this arrangement is that we save two multiplications per 1-D IDCT,
  119.  * because the y0 and y4 inputs need not be divided by sqrt(N).
  120.  *
  121.  * We have to do addition and subtraction of the integer inputs, which
  122.  * is no problem, and multiplication by fractional constants, which is
  123.  * a problem to do in integer arithmetic.  We multiply all the constants
  124.  * by CONST_SCALE and convert them to integer constants (thus retaining
  125.  * CONST_BITS bits of precision in the constants).  After doing a
  126.  * multiplication we have to divide the product by CONST_SCALE, with proper
  127.  * rounding, to produce the correct output.  This division can be done
  128.  * cheaply as a right shift of CONST_BITS bits.  We postpone shifting
  129.  * as long as possible so that partial sums can be added together with
  130.  * full fractional precision.
  131.  *
  132.  * The outputs of the first pass are scaled up by PASS1_BITS bits so that
  133.  * they are represented to better-than-integral precision.  These outputs
  134.  * require BITS_IN_JSAMPLE + PASS1_BITS + 3 bits; this fits in a 16-bit word
  135.  * with the recommended scaling.  (To scale up 12-bit sample data further, an
  136.  * intermediate INT32 array would be needed.)
  137.  *
  138.  * To avoid overflow of the 32-bit intermediate results in pass 2, we must
  139.  * have BITS_IN_JSAMPLE + CONST_BITS + PASS1_BITS <= 26.  Error analysis
  140.  * shows that the values given below are the most effective.
  141.  */
  142.  
  143. #ifdef EIGHT_BIT_SAMPLES
  144. #define PASS1_BITS  2
  145. #else
  146. #define PASS1_BITS  1        /* lose a little precision to avoid overflow */
  147. #endif
  148.  
  149. #define ONE    ((INT32) 1)
  150.  
  151. #define CONST_SCALE (ONE << CONST_BITS)
  152.  
  153. /* Convert a positive real constant to an integer scaled by CONST_SCALE.
  154.  * IMPORTANT: if your compiler doesn't do this arithmetic at compile time,
  155.  * you will pay a significant penalty in run time.  In that case, figure
  156.  * the correct integer constant values and insert them by hand.
  157.  */
  158.  
  159. #define FIX(x)    ((INT32) ((x) * CONST_SCALE + 0.5))
  160.  
  161. /* Descale and correctly round an INT32 value that's scaled by N bits.
  162.  * We assume RIGHT_SHIFT rounds towards minus infinity, so adding
  163.  * the fudge factor is correct for either sign of X.
  164.  */
  165.  
  166. #define DESCALE(x,n)  RIGHT_SHIFT((x) + (ONE << ((n)-1)), n)
  167.  
  168. /* Multiply an INT32 variable by an INT32 constant to yield an INT32 result.
  169.  * For 8-bit samples with the recommended scaling, all the variable
  170.  * and constant values involved are no more than 16 bits wide, so a
  171.  * 16x16->32 bit multiply can be used instead of a full 32x32 multiply;
  172.  * this provides a useful speedup on many machines.
  173.  * There is no way to specify a 16x16->32 multiply in portable C, but
  174.  * some C compilers will do the right thing if you provide the correct
  175.  * combination of casts.
  176.  * NB: for 12-bit samples, a full 32-bit multiplication will be needed.
  177.  */
  178.  
  179. #ifdef EIGHT_BIT_SAMPLES
  180. #ifdef SHORTxSHORT_32        /* may work if 'int' is 32 bits */
  181. #define MULTIPLY(var,const)  (((INT16) (var)) * ((INT16) (const)))
  182. #endif
  183. #ifdef SHORTxLCONST_32        /* known to work with Microsoft C 6.0 */
  184. #define MULTIPLY(var,const)  (((INT16) (var)) * ((INT32) (const)))
  185. #endif
  186. #endif
  187.  
  188. #ifndef MULTIPLY        /* default definition */
  189. #define MULTIPLY(var,const)  ((var) * (const))
  190. #endif
  191.  
  192. /* Precomputed idct value arrays. */
  193.  
  194. static DCTELEM PreIDCT[64][64];
  195.  
  196. /* Pre compute singleton coefficient IDCT values. */
  197. void
  198. init_pre_idct() {
  199.   int i;
  200.   void j_rev_dct();
  201.  
  202.   for (i=0; i<64; i++) {
  203.     memset((char *) PreIDCT[i], 0, 64*sizeof(DCTELEM));
  204.     PreIDCT[i][i] = 2048;
  205.     j_rev_dct(PreIDCT[i]);
  206.   }
  207. }
  208.  
  209. #ifndef ORIG_DCT
  210.   
  211.  
  212. /*
  213.  * Perform the inverse DCT on one block of coefficients.
  214.  */
  215.  
  216. void
  217. j_rev_dct_sparse (data, pos)
  218.      DCTBLOCK data;
  219.      int pos;
  220. {
  221.   register DCTELEM *dataptr;
  222.   short int val;
  223.   DCTELEM *ndataptr;
  224.   int scale, coeff, rr;
  225.   register int *dp;
  226.   register int v;
  227.  
  228.   /* If DC Coefficient. */
  229.   
  230.   if (pos == 0) {
  231.     dp = (int *)data;
  232.     v = *data;
  233.     /* Compute 32 bit value to assign.  This speeds things up a bit */
  234.     if (v < 0) val = (v-3)>>3;
  235.     else val = (v+4)>>3;
  236.     v = val | (val << 16);
  237.     dp[0] = v;      dp[1] = v;      dp[2] = v;      dp[3] = v;
  238.     dp[4] = v;      dp[5] = v;      dp[6] = v;      dp[7] = v;
  239.     dp[8] = v;      dp[9] = v;      dp[10] = v;      dp[11] = v;
  240.     dp[12] = v;      dp[13] = v;      dp[14] = v;      dp[15] = v;
  241.     dp[16] = v;      dp[17] = v;      dp[18] = v;      dp[19] = v;
  242.     dp[20] = v;      dp[21] = v;      dp[22] = v;      dp[23] = v;
  243.     dp[24] = v;      dp[25] = v;      dp[26] = v;      dp[27] = v;
  244.     dp[28] = v;      dp[29] = v;      dp[30] = v;      dp[31] = v;
  245.     return;
  246.   }
  247.   
  248.   /* Some other coefficient. */
  249.   dataptr = (DCTELEM *)data;
  250.   coeff = dataptr[pos];
  251.   ndataptr = PreIDCT[pos];
  252.  
  253.   for (rr=0; rr<4; rr++) {
  254.     dataptr[0] = (ndataptr[0] * coeff) >> (CONST_BITS-2);
  255.     dataptr[1] = (ndataptr[1] * coeff) >> (CONST_BITS-2);
  256.     dataptr[2] = (ndataptr[2] * coeff) >> (CONST_BITS-2);
  257.     dataptr[3] = (ndataptr[3] * coeff) >> (CONST_BITS-2);
  258.     dataptr[4] = (ndataptr[4] * coeff) >> (CONST_BITS-2);
  259.     dataptr[5] = (ndataptr[5] * coeff) >> (CONST_BITS-2);
  260.     dataptr[6] = (ndataptr[6] * coeff) >> (CONST_BITS-2);
  261.     dataptr[7] = (ndataptr[7] * coeff) >> (CONST_BITS-2);
  262.     dataptr[8] = (ndataptr[8] * coeff) >> (CONST_BITS-2);
  263.     dataptr[9] = (ndataptr[9] * coeff) >> (CONST_BITS-2);
  264.     dataptr[10] = (ndataptr[10] * coeff) >> (CONST_BITS-2);
  265.     dataptr[11] = (ndataptr[11] * coeff) >> (CONST_BITS-2);
  266.     dataptr[12] = (ndataptr[12] * coeff) >> (CONST_BITS-2);
  267.     dataptr[13] = (ndataptr[13] * coeff) >> (CONST_BITS-2);
  268.     dataptr[14] = (ndataptr[14] * coeff) >> (CONST_BITS-2);
  269.     dataptr[15] = (ndataptr[15] * coeff) >> (CONST_BITS-2);
  270.     dataptr += 16;
  271.     ndataptr += 16;
  272.   }
  273.   return;
  274. }
  275.  
  276.  
  277. void
  278. j_rev_dct (data)
  279.      DCTBLOCK data;
  280. {
  281.   INT32 tmp0, tmp1, tmp2, tmp3;
  282.   INT32 tmp10, tmp11, tmp12, tmp13;
  283.   INT32 z1, z2, z3, z4, z5;
  284.   INT32 d0, d1, d2, d3, d4, d5, d6, d7;
  285.   register DCTELEM *dataptr;
  286.   int rowctr;
  287.   SHIFT_TEMPS
  288.    
  289.   /* Pass 1: process rows. */
  290.   /* Note results are scaled up by sqrt(8) compared to a true IDCT; */
  291.   /* furthermore, we scale the results by 2**PASS1_BITS. */
  292.  
  293.   dataptr = data;
  294.  
  295.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  296.     /* Due to quantization, we will usually find that many of the input
  297.      * coefficients are zero, especially the AC terms.  We can exploit this
  298.      * by short-circuiting the IDCT calculation for any row in which all
  299.      * the AC terms are zero.  In that case each output is equal to the
  300.      * DC coefficient (with scale factor as needed).
  301.      * With typical images and quantization tables, half or more of the
  302.      * row DCT calculations can be simplified this way.
  303.      */
  304.  
  305.     register int *idataptr = (int*)dataptr;
  306.     d0 = dataptr[0];
  307.     d1 = dataptr[1];
  308.     if ((d1 == 0) && (idataptr[1] | idataptr[2] | idataptr[3]) == 0) {
  309.       /* AC terms all zero */
  310.       if (d0) {
  311.       /* Compute a 32 bit value to assign. */
  312.       DCTELEM dcval = (DCTELEM) (d0 << PASS1_BITS);
  313.       register int v = (dcval & 0xffff) | ((dcval << 16) & 0xffff0000);
  314.       
  315.       idataptr[0] = v;
  316.       idataptr[1] = v;
  317.       idataptr[2] = v;
  318.       idataptr[3] = v;
  319.       }
  320.       
  321.       dataptr += DCTSIZE;    /* advance pointer to next row */
  322.       continue;
  323.     }
  324.     d2 = dataptr[2];
  325.     d3 = dataptr[3];
  326.     d4 = dataptr[4];
  327.     d5 = dataptr[5];
  328.     d6 = dataptr[6];
  329.     d7 = dataptr[7];
  330.  
  331.     /* Even part: reverse the even part of the forward DCT. */
  332.     /* The rotator is sqrt(2)*c(-6). */
  333.     if (d6) {
  334.     if (d4) {
  335.         if (d2) {
  336.         if (d0) {
  337.             /* d0 != 0, d2 != 0, d4 != 0, d6 != 0 */
  338.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  339.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  340.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  341.  
  342.             tmp0 = (d0 + d4) << CONST_BITS;
  343.             tmp1 = (d0 - d4) << CONST_BITS;
  344.  
  345.             tmp10 = tmp0 + tmp3;
  346.             tmp13 = tmp0 - tmp3;
  347.             tmp11 = tmp1 + tmp2;
  348.             tmp12 = tmp1 - tmp2;
  349.         } else {
  350.             /* d0 == 0, d2 != 0, d4 != 0, d6 != 0 */
  351.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  352.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  353.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  354.  
  355.             tmp0 = d4 << CONST_BITS;
  356.  
  357.             tmp10 = tmp0 + tmp3;
  358.             tmp13 = tmp0 - tmp3;
  359.             tmp11 = tmp2 - tmp0;
  360.             tmp12 = -(tmp0 + tmp2);
  361.         }
  362.         } else {
  363.         if (d0) {
  364.             /* d0 != 0, d2 == 0, d4 != 0, d6 != 0 */
  365.             tmp2 = MULTIPLY(d6, - FIX(1.306562965));
  366.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  367.  
  368.             tmp0 = (d0 + d4) << CONST_BITS;
  369.             tmp1 = (d0 - d4) << CONST_BITS;
  370.  
  371.             tmp10 = tmp0 + tmp3;
  372.             tmp13 = tmp0 - tmp3;
  373.             tmp11 = tmp1 + tmp2;
  374.             tmp12 = tmp1 - tmp2;
  375.         } else {
  376.             /* d0 == 0, d2 == 0, d4 != 0, d6 != 0 */
  377.             tmp2 = MULTIPLY(d6, -FIX(1.306562965));
  378.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  379.  
  380.             tmp0 = d4 << CONST_BITS;
  381.  
  382.             tmp10 = tmp0 + tmp3;
  383.             tmp13 = tmp0 - tmp3;
  384.             tmp11 = tmp2 - tmp0;
  385.             tmp12 = -(tmp0 + tmp2);
  386.         }
  387.         }
  388.     } else {
  389.         if (d2) {
  390.         if (d0) {
  391.             /* d0 != 0, d2 != 0, d4 == 0, d6 != 0 */
  392.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  393.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  394.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  395.  
  396.             tmp0 = d0 << CONST_BITS;
  397.  
  398.             tmp10 = tmp0 + tmp3;
  399.             tmp13 = tmp0 - tmp3;
  400.             tmp11 = tmp0 + tmp2;
  401.             tmp12 = tmp0 - tmp2;
  402.         } else {
  403.             /* d0 == 0, d2 != 0, d4 == 0, d6 != 0 */
  404.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  405.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  406.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  407.  
  408.             tmp10 = tmp3;
  409.             tmp13 = -tmp3;
  410.             tmp11 = tmp2;
  411.             tmp12 = -tmp2;
  412.         }
  413.         } else {
  414.         if (d0) {
  415.             /* d0 != 0, d2 == 0, d4 == 0, d6 != 0 */
  416.             tmp2 = MULTIPLY(d6, - FIX(1.306562965));
  417.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  418.  
  419.             tmp0 = d0 << CONST_BITS;
  420.  
  421.             tmp10 = tmp0 + tmp3;
  422.             tmp13 = tmp0 - tmp3;
  423.             tmp11 = tmp0 + tmp2;
  424.             tmp12 = tmp0 - tmp2;
  425.         } else {
  426.             /* d0 == 0, d2 == 0, d4 == 0, d6 != 0 */
  427.             tmp2 = MULTIPLY(d6, - FIX(1.306562965));
  428.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  429.  
  430.             tmp10 = tmp3;
  431.             tmp13 = -tmp3;
  432.             tmp11 = tmp2;
  433.             tmp12 = -tmp2;
  434.         }
  435.         }
  436.     }
  437.     } else {
  438.     if (d4) {
  439.         if (d2) {
  440.         if (d0) {
  441.             /* d0 != 0, d2 != 0, d4 != 0, d6 == 0 */
  442.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  443.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  444.  
  445.             tmp0 = (d0 + d4) << CONST_BITS;
  446.             tmp1 = (d0 - d4) << CONST_BITS;
  447.  
  448.             tmp10 = tmp0 + tmp3;
  449.             tmp13 = tmp0 - tmp3;
  450.             tmp11 = tmp1 + tmp2;
  451.             tmp12 = tmp1 - tmp2;
  452.         } else {
  453.             /* d0 == 0, d2 != 0, d4 != 0, d6 == 0 */
  454.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  455.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  456.  
  457.             tmp0 = d4 << CONST_BITS;
  458.  
  459.             tmp10 = tmp0 + tmp3;
  460.             tmp13 = tmp0 - tmp3;
  461.             tmp11 = tmp2 - tmp0;
  462.             tmp12 = -(tmp0 + tmp2);
  463.         }
  464.         } else {
  465.         if (d0) {
  466.             /* d0 != 0, d2 == 0, d4 != 0, d6 == 0 */
  467.             tmp10 = tmp13 = (d0 + d4) << CONST_BITS;
  468.             tmp11 = tmp12 = (d0 - d4) << CONST_BITS;
  469.         } else {
  470.             /* d0 == 0, d2 == 0, d4 != 0, d6 == 0 */
  471.             tmp10 = tmp13 = d4 << CONST_BITS;
  472.             tmp11 = tmp12 = -tmp10;
  473.         }
  474.         }
  475.     } else {
  476.         if (d2) {
  477.         if (d0) {
  478.             /* d0 != 0, d2 != 0, d4 == 0, d6 == 0 */
  479.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  480.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  481.  
  482.             tmp0 = d0 << CONST_BITS;
  483.  
  484.             tmp10 = tmp0 + tmp3;
  485.             tmp13 = tmp0 - tmp3;
  486.             tmp11 = tmp0 + tmp2;
  487.             tmp12 = tmp0 - tmp2;
  488.         } else {
  489.             /* d0 == 0, d2 != 0, d4 == 0, d6 == 0 */
  490.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  491.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  492.  
  493.             tmp10 = tmp3;
  494.             tmp13 = -tmp3;
  495.             tmp11 = tmp2;
  496.             tmp12 = -tmp2;
  497.         }
  498.         } else {
  499.         if (d0) {
  500.             /* d0 != 0, d2 == 0, d4 == 0, d6 == 0 */
  501.             tmp10 = tmp13 = tmp11 = tmp12 = d0 << CONST_BITS;
  502.         } else {
  503.             /* d0 == 0, d2 == 0, d4 == 0, d6 == 0 */
  504.             tmp10 = tmp13 = tmp11 = tmp12 = 0;
  505.         }
  506.         }
  507.     }
  508.     }
  509.  
  510.  
  511.     /* Odd part per figure 8; the matrix is unitary and hence its
  512.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  513.      */
  514.  
  515.     if (d7) {
  516.     if (d5) {
  517.         if (d3) {
  518.         if (d1) {
  519.             /* d1 != 0, d3 != 0, d5 != 0, d7 != 0 */
  520.             z1 = d7 + d1;
  521.             z2 = d5 + d3;
  522.             z3 = d7 + d3;
  523.             z4 = d5 + d1;
  524.             z5 = MULTIPLY(z3 + z4, FIX(1.175875602));
  525.             
  526.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  527.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  528.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  529.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  530.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  531.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  532.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  533.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  534.             
  535.             z3 += z5;
  536.             z4 += z5;
  537.             
  538.             tmp0 += z1 + z3;
  539.             tmp1 += z2 + z4;
  540.             tmp2 += z2 + z3;
  541.             tmp3 += z1 + z4;
  542.         } else {
  543.             /* d1 == 0, d3 != 0, d5 != 0, d7 != 0 */
  544.             z1 = d7;
  545.             z2 = d5 + d3;
  546.             z3 = d7 + d3;
  547.             z5 = MULTIPLY(z3 + d5, FIX(1.175875602));
  548.             
  549.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  550.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  551.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  552.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  553.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  554.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  555.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  556.             
  557.             z3 += z5;
  558.             z4 += z5;
  559.             
  560.             tmp0 += z1 + z3;
  561.             tmp1 += z2 + z4;
  562.             tmp2 += z2 + z3;
  563.             tmp3 = z1 + z4;
  564.         }
  565.         } else {
  566.         if (d1) {
  567.             /* d1 != 0, d3 == 0, d5 != 0, d7 != 0 */
  568.             z1 = d7 + d1;
  569.             z2 = d5;
  570.             z3 = d7;
  571.             z4 = d5 + d1;
  572.             z5 = MULTIPLY(z3 + z4, FIX(1.175875602));
  573.             
  574.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  575.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  576.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  577.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  578.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  579.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  580.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  581.             
  582.             z3 += z5;
  583.             z4 += z5;
  584.             
  585.             tmp0 += z1 + z3;
  586.             tmp1 += z2 + z4;
  587.             tmp2 = z2 + z3;
  588.             tmp3 += z1 + z4;
  589.         } else {
  590.             /* d1 == 0, d3 == 0, d5 != 0, d7 != 0 */
  591.             tmp0 = MULTIPLY(d7, - FIX(0.601344887)); 
  592.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  593.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  594.             tmp1 = MULTIPLY(d5, - FIX(0.509795578));
  595.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  596.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  597.             z5 = MULTIPLY(d5 + d7, FIX(1.175875602));
  598.             
  599.             z3 += z5;
  600.             z4 += z5;
  601.             
  602.             tmp0 += z3;
  603.             tmp1 += z4;
  604.             tmp2 = z2 + z3;
  605.             tmp3 = z1 + z4;
  606.         }
  607.         }
  608.     } else {
  609.         if (d3) {
  610.         if (d1) {
  611.             /* d1 != 0, d3 != 0, d5 == 0, d7 != 0 */
  612.             z1 = d7 + d1;
  613.             z3 = d7 + d3;
  614.             z5 = MULTIPLY(z3 + d1, FIX(1.175875602));
  615.             
  616.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  617.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  618.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  619.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  620.             z2 = MULTIPLY(d3, - FIX(2.562915447));
  621.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  622.             z4 = MULTIPLY(d1, - FIX(0.390180644));
  623.             
  624.             z3 += z5;
  625.             z4 += z5;
  626.             
  627.             tmp0 += z1 + z3;
  628.             tmp1 = z2 + z4;
  629.             tmp2 += z2 + z3;
  630.             tmp3 += z1 + z4;
  631.         } else {
  632.             /* d1 == 0, d3 != 0, d5 == 0, d7 != 0 */
  633.             z3 = d7 + d3;
  634.             
  635.             tmp0 = MULTIPLY(d7, - FIX(0.601344887)); 
  636.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  637.             tmp2 = MULTIPLY(d3, FIX(0.509795579));
  638.             z2 = MULTIPLY(d3, - FIX(2.562915447));
  639.             z5 = MULTIPLY(z3, FIX(1.175875602));
  640.             z3 = MULTIPLY(z3, - FIX(0.785694958));
  641.             
  642.             tmp0 += z3;
  643.             tmp1 = z2 + z5;
  644.             tmp2 += z3;
  645.             tmp3 = z1 + z5;
  646.         }
  647.         } else {
  648.         if (d1) {
  649.             /* d1 != 0, d3 == 0, d5 == 0, d7 != 0 */
  650.             z1 = d7 + d1;
  651.             z5 = MULTIPLY(z1, FIX(1.175875602));
  652.  
  653.             z1 = MULTIPLY(z1, FIX(0.275899379));
  654.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  655.             tmp0 = MULTIPLY(d7, - FIX(1.662939224)); 
  656.             z4 = MULTIPLY(d1, - FIX(0.390180644));
  657.             tmp3 = MULTIPLY(d1, FIX(1.111140466));
  658.  
  659.             tmp0 += z1;
  660.             tmp1 = z4 + z5;
  661.             tmp2 = z3 + z5;
  662.             tmp3 += z1;
  663.         } else {
  664.             /* d1 == 0, d3 == 0, d5 == 0, d7 != 0 */
  665.             tmp0 = MULTIPLY(d7, - FIX(1.387039845));
  666.             tmp1 = MULTIPLY(d7, FIX(1.175875602));
  667.             tmp2 = MULTIPLY(d7, - FIX(0.785694958));
  668.             tmp3 = MULTIPLY(d7, FIX(0.275899379));
  669.         }
  670.         }
  671.     }
  672.     } else {
  673.     if (d5) {
  674.         if (d3) {
  675.         if (d1) {
  676.             /* d1 != 0, d3 != 0, d5 != 0, d7 == 0 */
  677.             z2 = d5 + d3;
  678.             z4 = d5 + d1;
  679.             z5 = MULTIPLY(d3 + z4, FIX(1.175875602));
  680.             
  681.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  682.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  683.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  684.             z1 = MULTIPLY(d1, - FIX(0.899976223));
  685.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  686.             z3 = MULTIPLY(d3, - FIX(1.961570560));
  687.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  688.             
  689.             z3 += z5;
  690.             z4 += z5;
  691.             
  692.             tmp0 = z1 + z3;
  693.             tmp1 += z2 + z4;
  694.             tmp2 += z2 + z3;
  695.             tmp3 += z1 + z4;
  696.         } else {
  697.             /* d1 == 0, d3 != 0, d5 != 0, d7 == 0 */
  698.             z2 = d5 + d3;
  699.             
  700.             z5 = MULTIPLY(z2, FIX(1.175875602));
  701.             tmp1 = MULTIPLY(d5, FIX(1.662939225));
  702.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  703.             z2 = MULTIPLY(z2, - FIX(1.387039845));
  704.             tmp2 = MULTIPLY(d3, FIX(1.111140466));
  705.             z3 = MULTIPLY(d3, - FIX(1.961570560));
  706.             
  707.             tmp0 = z3 + z5;
  708.             tmp1 += z2;
  709.             tmp2 += z2;
  710.             tmp3 = z4 + z5;
  711.         }
  712.         } else {
  713.         if (d1) {
  714.             /* d1 != 0, d3 == 0, d5 != 0, d7 == 0 */
  715.             z4 = d5 + d1;
  716.             
  717.             z5 = MULTIPLY(z4, FIX(1.175875602));
  718.             z1 = MULTIPLY(d1, - FIX(0.899976223));
  719.             tmp3 = MULTIPLY(d1, FIX(0.601344887));
  720.             tmp1 = MULTIPLY(d5, - FIX(0.509795578));
  721.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  722.             z4 = MULTIPLY(z4, FIX(0.785694958));
  723.             
  724.             tmp0 = z1 + z5;
  725.             tmp1 += z4;
  726.             tmp2 = z2 + z5;
  727.             tmp3 += z4;
  728.         } else {
  729.             /* d1 == 0, d3 == 0, d5 != 0, d7 == 0 */
  730.             tmp0 = MULTIPLY(d5, FIX(1.175875602));
  731.             tmp1 = MULTIPLY(d5, FIX(0.275899380));
  732.             tmp2 = MULTIPLY(d5, - FIX(1.387039845));
  733.             tmp3 = MULTIPLY(d5, FIX(0.785694958));
  734.         }
  735.         }
  736.     } else {
  737.         if (d3) {
  738.         if (d1) {
  739.             /* d1 != 0, d3 != 0, d5 == 0, d7 == 0 */
  740.             z5 = d1 + d3;
  741.             tmp3 = MULTIPLY(d1, FIX(0.211164243));
  742.             tmp2 = MULTIPLY(d3, - FIX(1.451774981));
  743.             z1 = MULTIPLY(d1, FIX(1.061594337));
  744.             z2 = MULTIPLY(d3, - FIX(2.172734803));
  745.             z4 = MULTIPLY(z5, FIX(0.785694958));
  746.             z5 = MULTIPLY(z5, FIX(1.175875602));
  747.             
  748.             tmp0 = z1 - z4;
  749.             tmp1 = z2 + z4;
  750.             tmp2 += z5;
  751.             tmp3 += z5;
  752.         } else {
  753.             /* d1 == 0, d3 != 0, d5 == 0, d7 == 0 */
  754.             tmp0 = MULTIPLY(d3, - FIX(0.785694958));
  755.             tmp1 = MULTIPLY(d3, - FIX(1.387039845));
  756.             tmp2 = MULTIPLY(d3, - FIX(0.275899379));
  757.             tmp3 = MULTIPLY(d3, FIX(1.175875602));
  758.         }
  759.         } else {
  760.         if (d1) {
  761.             /* d1 != 0, d3 == 0, d5 == 0, d7 == 0 */
  762.             tmp0 = MULTIPLY(d1, FIX(0.275899379));
  763.             tmp1 = MULTIPLY(d1, FIX(0.785694958));
  764.             tmp2 = MULTIPLY(d1, FIX(1.175875602));
  765.             tmp3 = MULTIPLY(d1, FIX(1.387039845));
  766.         } else {
  767.             /* d1 == 0, d3 == 0, d5 == 0, d7 == 0 */
  768.             tmp0 = tmp1 = tmp2 = tmp3 = 0;
  769.         }
  770.         }
  771.     }
  772.     }
  773.  
  774.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  775.  
  776.     dataptr[0] = (DCTELEM) DESCALE(tmp10 + tmp3, CONST_BITS-PASS1_BITS);
  777.     dataptr[7] = (DCTELEM) DESCALE(tmp10 - tmp3, CONST_BITS-PASS1_BITS);
  778.     dataptr[1] = (DCTELEM) DESCALE(tmp11 + tmp2, CONST_BITS-PASS1_BITS);
  779.     dataptr[6] = (DCTELEM) DESCALE(tmp11 - tmp2, CONST_BITS-PASS1_BITS);
  780.     dataptr[2] = (DCTELEM) DESCALE(tmp12 + tmp1, CONST_BITS-PASS1_BITS);
  781.     dataptr[5] = (DCTELEM) DESCALE(tmp12 - tmp1, CONST_BITS-PASS1_BITS);
  782.     dataptr[3] = (DCTELEM) DESCALE(tmp13 + tmp0, CONST_BITS-PASS1_BITS);
  783.     dataptr[4] = (DCTELEM) DESCALE(tmp13 - tmp0, CONST_BITS-PASS1_BITS);
  784.  
  785.     dataptr += DCTSIZE;        /* advance pointer to next row */
  786.   }
  787.  
  788.   /* Pass 2: process columns. */
  789.   /* Note that we must descale the results by a factor of 8 == 2**3, */
  790.   /* and also undo the PASS1_BITS scaling. */
  791.  
  792.   dataptr = data;
  793.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  794.     /* Columns of zeroes can be exploited in the same way as we did with rows.
  795.      * However, the row calculation has created many nonzero AC terms, so the
  796.      * simplification applies less often (typically 5% to 10% of the time).
  797.      * On machines with very fast multiplication, it's possible that the
  798.      * test takes more time than it's worth.  In that case this section
  799.      * may be commented out.
  800.      */
  801.  
  802.     d0 = dataptr[DCTSIZE*0];
  803.     d1 = dataptr[DCTSIZE*1];
  804.     d2 = dataptr[DCTSIZE*2];
  805.     d3 = dataptr[DCTSIZE*3];
  806.     d4 = dataptr[DCTSIZE*4];
  807.     d5 = dataptr[DCTSIZE*5];
  808.     d6 = dataptr[DCTSIZE*6];
  809.     d7 = dataptr[DCTSIZE*7];
  810.  
  811.     /* Even part: reverse the even part of the forward DCT. */
  812.     /* The rotator is sqrt(2)*c(-6). */
  813.     if (d6) {
  814.     if (d4) {
  815.         if (d2) {
  816.         if (d0) {
  817.             /* d0 != 0, d2 != 0, d4 != 0, d6 != 0 */
  818.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  819.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  820.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  821.  
  822.             tmp0 = (d0 + d4) << CONST_BITS;
  823.             tmp1 = (d0 - d4) << CONST_BITS;
  824.  
  825.             tmp10 = tmp0 + tmp3;
  826.             tmp13 = tmp0 - tmp3;
  827.             tmp11 = tmp1 + tmp2;
  828.             tmp12 = tmp1 - tmp2;
  829.         } else {
  830.             /* d0 == 0, d2 != 0, d4 != 0, d6 != 0 */
  831.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  832.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  833.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  834.  
  835.             tmp0 = d4 << CONST_BITS;
  836.  
  837.             tmp10 = tmp0 + tmp3;
  838.             tmp13 = tmp0 - tmp3;
  839.             tmp11 = tmp2 - tmp0;
  840.             tmp12 = -(tmp0 + tmp2);
  841.         }
  842.         } else {
  843.         if (d0) {
  844.             /* d0 != 0, d2 == 0, d4 != 0, d6 != 0 */
  845.             tmp2 = MULTIPLY(d6, - FIX(1.306562965));
  846.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  847.  
  848.             tmp0 = (d0 + d4) << CONST_BITS;
  849.             tmp1 = (d0 - d4) << CONST_BITS;
  850.  
  851.             tmp10 = tmp0 + tmp3;
  852.             tmp13 = tmp0 - tmp3;
  853.             tmp11 = tmp1 + tmp2;
  854.             tmp12 = tmp1 - tmp2;
  855.         } else {
  856.             /* d0 == 0, d2 == 0, d4 != 0, d6 != 0 */
  857.             tmp2 = MULTIPLY(d6, -FIX(1.306562965));
  858.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  859.  
  860.             tmp0 = d4 << CONST_BITS;
  861.  
  862.             tmp10 = tmp0 + tmp3;
  863.             tmp13 = tmp0 - tmp3;
  864.             tmp11 = tmp2 - tmp0;
  865.             tmp12 = -(tmp0 + tmp2);
  866.         }
  867.         }
  868.     } else {
  869.         if (d2) {
  870.         if (d0) {
  871.             /* d0 != 0, d2 != 0, d4 == 0, d6 != 0 */
  872.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  873.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  874.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  875.  
  876.             tmp0 = d0 << CONST_BITS;
  877.  
  878.             tmp10 = tmp0 + tmp3;
  879.             tmp13 = tmp0 - tmp3;
  880.             tmp11 = tmp0 + tmp2;
  881.             tmp12 = tmp0 - tmp2;
  882.         } else {
  883.             /* d0 == 0, d2 != 0, d4 == 0, d6 != 0 */
  884.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  885.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  886.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  887.  
  888.             tmp10 = tmp3;
  889.             tmp13 = -tmp3;
  890.             tmp11 = tmp2;
  891.             tmp12 = -tmp2;
  892.         }
  893.         } else {
  894.         if (d0) {
  895.             /* d0 != 0, d2 == 0, d4 == 0, d6 != 0 */
  896.             tmp2 = MULTIPLY(d6, - FIX(1.306562965));
  897.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  898.  
  899.             tmp0 = d0 << CONST_BITS;
  900.  
  901.             tmp10 = tmp0 + tmp3;
  902.             tmp13 = tmp0 - tmp3;
  903.             tmp11 = tmp0 + tmp2;
  904.             tmp12 = tmp0 - tmp2;
  905.         } else {
  906.             /* d0 == 0, d2 == 0, d4 == 0, d6 != 0 */
  907.             tmp2 = MULTIPLY(d6, - FIX(1.306562965));
  908.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  909.  
  910.             tmp10 = tmp3;
  911.             tmp13 = -tmp3;
  912.             tmp11 = tmp2;
  913.             tmp12 = -tmp2;
  914.         }
  915.         }
  916.     }
  917.     } else {
  918.     if (d4) {
  919.         if (d2) {
  920.         if (d0) {
  921.             /* d0 != 0, d2 != 0, d4 != 0, d6 == 0 */
  922.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  923.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  924.  
  925.             tmp0 = (d0 + d4) << CONST_BITS;
  926.             tmp1 = (d0 - d4) << CONST_BITS;
  927.  
  928.             tmp10 = tmp0 + tmp3;
  929.             tmp13 = tmp0 - tmp3;
  930.             tmp11 = tmp1 + tmp2;
  931.             tmp12 = tmp1 - tmp2;
  932.         } else {
  933.             /* d0 == 0, d2 != 0, d4 != 0, d6 == 0 */
  934.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  935.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  936.  
  937.             tmp0 = d4 << CONST_BITS;
  938.  
  939.             tmp10 = tmp0 + tmp3;
  940.             tmp13 = tmp0 - tmp3;
  941.             tmp11 = tmp2 - tmp0;
  942.             tmp12 = -(tmp0 + tmp2);
  943.         }
  944.         } else {
  945.         if (d0) {
  946.             /* d0 != 0, d2 == 0, d4 != 0, d6 == 0 */
  947.             tmp10 = tmp13 = (d0 + d4) << CONST_BITS;
  948.             tmp11 = tmp12 = (d0 - d4) << CONST_BITS;
  949.         } else {
  950.             /* d0 == 0, d2 == 0, d4 != 0, d6 == 0 */
  951.             tmp10 = tmp13 = d4 << CONST_BITS;
  952.             tmp11 = tmp12 = -tmp10;
  953.         }
  954.         }
  955.     } else {
  956.         if (d2) {
  957.         if (d0) {
  958.             /* d0 != 0, d2 != 0, d4 == 0, d6 == 0 */
  959.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  960.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  961.  
  962.             tmp0 = d0 << CONST_BITS;
  963.  
  964.             tmp10 = tmp0 + tmp3;
  965.             tmp13 = tmp0 - tmp3;
  966.             tmp11 = tmp0 + tmp2;
  967.             tmp12 = tmp0 - tmp2;
  968.         } else {
  969.             /* d0 == 0, d2 != 0, d4 == 0, d6 == 0 */
  970.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  971.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  972.  
  973.             tmp10 = tmp3;
  974.             tmp13 = -tmp3;
  975.             tmp11 = tmp2;
  976.             tmp12 = -tmp2;
  977.         }
  978.         } else {
  979.         if (d0) {
  980.             /* d0 != 0, d2 == 0, d4 == 0, d6 == 0 */
  981.             tmp10 = tmp13 = tmp11 = tmp12 = d0 << CONST_BITS;
  982.         } else {
  983.             /* d0 == 0, d2 == 0, d4 == 0, d6 == 0 */
  984.             tmp10 = tmp13 = tmp11 = tmp12 = 0;
  985.         }
  986.         }
  987.     }
  988.     }
  989.  
  990.     /* Odd part per figure 8; the matrix is unitary and hence its
  991.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  992.      */
  993.     if (d7) {
  994.     if (d5) {
  995.         if (d3) {
  996.         if (d1) {
  997.             /* d1 != 0, d3 != 0, d5 != 0, d7 != 0 */
  998.             z1 = d7 + d1;
  999.             z2 = d5 + d3;
  1000.             z3 = d7 + d3;
  1001.             z4 = d5 + d1;
  1002.             z5 = MULTIPLY(z3 + z4, FIX(1.175875602));
  1003.             
  1004.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  1005.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  1006.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  1007.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  1008.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  1009.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  1010.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  1011.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  1012.             
  1013.             z3 += z5;
  1014.             z4 += z5;
  1015.             
  1016.             tmp0 += z1 + z3;
  1017.             tmp1 += z2 + z4;
  1018.             tmp2 += z2 + z3;
  1019.             tmp3 += z1 + z4;
  1020.         } else {
  1021.             /* d1 == 0, d3 != 0, d5 != 0, d7 != 0 */
  1022.             z1 = d7;
  1023.             z2 = d5 + d3;
  1024.             z3 = d7 + d3;
  1025.             z5 = MULTIPLY(z3 + d5, FIX(1.175875602));
  1026.             
  1027.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  1028.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  1029.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  1030.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  1031.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  1032.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  1033.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  1034.             
  1035.             z3 += z5;
  1036.             z4 += z5;
  1037.             
  1038.             tmp0 += z1 + z3;
  1039.             tmp1 += z2 + z4;
  1040.             tmp2 += z2 + z3;
  1041.             tmp3 = z1 + z4;
  1042.         }
  1043.         } else {
  1044.         if (d1) {
  1045.             /* d1 != 0, d3 == 0, d5 != 0, d7 != 0 */
  1046.             z1 = d7 + d1;
  1047.             z2 = d5;
  1048.             z3 = d7;
  1049.             z4 = d5 + d1;
  1050.             z5 = MULTIPLY(z3 + z4, FIX(1.175875602));
  1051.             
  1052.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  1053.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  1054.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  1055.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  1056.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  1057.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  1058.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  1059.             
  1060.             z3 += z5;
  1061.             z4 += z5;
  1062.             
  1063.             tmp0 += z1 + z3;
  1064.             tmp1 += z2 + z4;
  1065.             tmp2 = z2 + z3;
  1066.             tmp3 += z1 + z4;
  1067.         } else {
  1068.             /* d1 == 0, d3 == 0, d5 != 0, d7 != 0 */
  1069.             tmp0 = MULTIPLY(d7, - FIX(0.601344887)); 
  1070.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  1071.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  1072.             tmp1 = MULTIPLY(d5, - FIX(0.509795578));
  1073.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  1074.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  1075.             z5 = MULTIPLY(d5 + d7, FIX(1.175875602));
  1076.             
  1077.             z3 += z5;
  1078.             z4 += z5;
  1079.             
  1080.             tmp0 += z3;
  1081.             tmp1 += z4;
  1082.             tmp2 = z2 + z3;
  1083.             tmp3 = z1 + z4;
  1084.         }
  1085.         }
  1086.     } else {
  1087.         if (d3) {
  1088.         if (d1) {
  1089.             /* d1 != 0, d3 != 0, d5 == 0, d7 != 0 */
  1090.             z1 = d7 + d1;
  1091.             z3 = d7 + d3;
  1092.             z5 = MULTIPLY(z3 + d1, FIX(1.175875602));
  1093.             
  1094.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  1095.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  1096.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  1097.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  1098.             z2 = MULTIPLY(d3, - FIX(2.562915447));
  1099.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  1100.             z4 = MULTIPLY(d1, - FIX(0.390180644));
  1101.             
  1102.             z3 += z5;
  1103.             z4 += z5;
  1104.             
  1105.             tmp0 += z1 + z3;
  1106.             tmp1 = z2 + z4;
  1107.             tmp2 += z2 + z3;
  1108.             tmp3 += z1 + z4;
  1109.         } else {
  1110.             /* d1 == 0, d3 != 0, d5 == 0, d7 != 0 */
  1111.             z3 = d7 + d3;
  1112.             
  1113.             tmp0 = MULTIPLY(d7, - FIX(0.601344887)); 
  1114.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  1115.             tmp2 = MULTIPLY(d3, FIX(0.509795579));
  1116.             z2 = MULTIPLY(d3, - FIX(2.562915447));
  1117.             z5 = MULTIPLY(z3, FIX(1.175875602));
  1118.             z3 = MULTIPLY(z3, - FIX(0.785694958));
  1119.             
  1120.             tmp0 += z3;
  1121.             tmp1 = z2 + z5;
  1122.             tmp2 += z3;
  1123.             tmp3 = z1 + z5;
  1124.         }
  1125.         } else {
  1126.         if (d1) {
  1127.             /* d1 != 0, d3 == 0, d5 == 0, d7 != 0 */
  1128.             z1 = d7 + d1;
  1129.             z5 = MULTIPLY(z1, FIX(1.175875602));
  1130.  
  1131.             z1 = MULTIPLY(z1, FIX(0.275899379));
  1132.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  1133.             tmp0 = MULTIPLY(d7, - FIX(1.662939224)); 
  1134.             z4 = MULTIPLY(d1, - FIX(0.390180644));
  1135.             tmp3 = MULTIPLY(d1, FIX(1.111140466));
  1136.  
  1137.             tmp0 += z1;
  1138.             tmp1 = z4 + z5;
  1139.             tmp2 = z3 + z5;
  1140.             tmp3 += z1;
  1141.         } else {
  1142.             /* d1 == 0, d3 == 0, d5 == 0, d7 != 0 */
  1143.             tmp0 = MULTIPLY(d7, - FIX(1.387039845));
  1144.             tmp1 = MULTIPLY(d7, FIX(1.175875602));
  1145.             tmp2 = MULTIPLY(d7, - FIX(0.785694958));
  1146.             tmp3 = MULTIPLY(d7, FIX(0.275899379));
  1147.         }
  1148.         }
  1149.     }
  1150.     } else {
  1151.     if (d5) {
  1152.         if (d3) {
  1153.         if (d1) {
  1154.             /* d1 != 0, d3 != 0, d5 != 0, d7 == 0 */
  1155.             z2 = d5 + d3;
  1156.             z4 = d5 + d1;
  1157.             z5 = MULTIPLY(d3 + z4, FIX(1.175875602));
  1158.             
  1159.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  1160.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  1161.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  1162.             z1 = MULTIPLY(d1, - FIX(0.899976223));
  1163.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  1164.             z3 = MULTIPLY(d3, - FIX(1.961570560));
  1165.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  1166.             
  1167.             z3 += z5;
  1168.             z4 += z5;
  1169.             
  1170.             tmp0 = z1 + z3;
  1171.             tmp1 += z2 + z4;
  1172.             tmp2 += z2 + z3;
  1173.             tmp3 += z1 + z4;
  1174.         } else {
  1175.             /* d1 == 0, d3 != 0, d5 != 0, d7 == 0 */
  1176.             z2 = d5 + d3;
  1177.             
  1178.             z5 = MULTIPLY(z2, FIX(1.175875602));
  1179.             tmp1 = MULTIPLY(d5, FIX(1.662939225));
  1180.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  1181.             z2 = MULTIPLY(z2, - FIX(1.387039845));
  1182.             tmp2 = MULTIPLY(d3, FIX(1.111140466));
  1183.             z3 = MULTIPLY(d3, - FIX(1.961570560));
  1184.             
  1185.             tmp0 = z3 + z5;
  1186.             tmp1 += z2;
  1187.             tmp2 += z2;
  1188.             tmp3 = z4 + z5;
  1189.         }
  1190.         } else {
  1191.         if (d1) {
  1192.             /* d1 != 0, d3 == 0, d5 != 0, d7 == 0 */
  1193.             z4 = d5 + d1;
  1194.             
  1195.             z5 = MULTIPLY(z4, FIX(1.175875602));
  1196.             z1 = MULTIPLY(d1, - FIX(0.899976223));
  1197.             tmp3 = MULTIPLY(d1, FIX(0.601344887));
  1198.             tmp1 = MULTIPLY(d5, - FIX(0.509795578));
  1199.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  1200.             z4 = MULTIPLY(z4, FIX(0.785694958));
  1201.             
  1202.             tmp0 = z1 + z5;
  1203.             tmp1 += z4;
  1204.             tmp2 = z2 + z5;
  1205.             tmp3 += z4;
  1206.         } else {
  1207.             /* d1 == 0, d3 == 0, d5 != 0, d7 == 0 */
  1208.             tmp0 = MULTIPLY(d5, FIX(1.175875602));
  1209.             tmp1 = MULTIPLY(d5, FIX(0.275899380));
  1210.             tmp2 = MULTIPLY(d5, - FIX(1.387039845));
  1211.             tmp3 = MULTIPLY(d5, FIX(0.785694958));
  1212.         }
  1213.         }
  1214.     } else {
  1215.         if (d3) {
  1216.         if (d1) {
  1217.             /* d1 != 0, d3 != 0, d5 == 0, d7 == 0 */
  1218.             z5 = d1 + d3;
  1219.             tmp3 = MULTIPLY(d1, FIX(0.211164243));
  1220.             tmp2 = MULTIPLY(d3, - FIX(1.451774981));
  1221.             z1 = MULTIPLY(d1, FIX(1.061594337));
  1222.             z2 = MULTIPLY(d3, - FIX(2.172734803));
  1223.             z4 = MULTIPLY(z5, FIX(0.785694958));
  1224.             z5 = MULTIPLY(z5, FIX(1.175875602));
  1225.             
  1226.             tmp0 = z1 - z4;
  1227.             tmp1 = z2 + z4;
  1228.             tmp2 += z5;
  1229.             tmp3 += z5;
  1230.         } else {
  1231.             /* d1 == 0, d3 != 0, d5 == 0, d7 == 0 */
  1232.             tmp0 = MULTIPLY(d3, - FIX(0.785694958));
  1233.             tmp1 = MULTIPLY(d3, - FIX(1.387039845));
  1234.             tmp2 = MULTIPLY(d3, - FIX(0.275899379));
  1235.             tmp3 = MULTIPLY(d3, FIX(1.175875602));
  1236.         }
  1237.         } else {
  1238.         if (d1) {
  1239.             /* d1 != 0, d3 == 0, d5 == 0, d7 == 0 */
  1240.             tmp0 = MULTIPLY(d1, FIX(0.275899379));
  1241.             tmp1 = MULTIPLY(d1, FIX(0.785694958));
  1242.             tmp2 = MULTIPLY(d1, FIX(1.175875602));
  1243.             tmp3 = MULTIPLY(d1, FIX(1.387039845));
  1244.         } else {
  1245.             /* d1 == 0, d3 == 0, d5 == 0, d7 == 0 */
  1246.             tmp0 = tmp1 = tmp2 = tmp3 = 0;
  1247.         }
  1248.         }
  1249.     }
  1250.     }
  1251.  
  1252.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  1253.  
  1254.     dataptr[DCTSIZE*0] = (DCTELEM) DESCALE(tmp10 + tmp3,
  1255.                        CONST_BITS+PASS1_BITS+3);
  1256.     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp10 - tmp3,
  1257.                        CONST_BITS+PASS1_BITS+3);
  1258.     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp11 + tmp2,
  1259.                        CONST_BITS+PASS1_BITS+3);
  1260.     dataptr[DCTSIZE*6] = (DCTELEM) DESCALE(tmp11 - tmp2,
  1261.                        CONST_BITS+PASS1_BITS+3);
  1262.     dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(tmp12 + tmp1,
  1263.                        CONST_BITS+PASS1_BITS+3);
  1264.     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp12 - tmp1,
  1265.                        CONST_BITS+PASS1_BITS+3);
  1266.     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp13 + tmp0,
  1267.                        CONST_BITS+PASS1_BITS+3);
  1268.     dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(tmp13 - tmp0,
  1269.                        CONST_BITS+PASS1_BITS+3);
  1270.     
  1271.     dataptr++;            /* advance pointer to next column */
  1272.   }
  1273. }
  1274.  
  1275. #else
  1276.  
  1277.  
  1278. void
  1279. j_rev_dct_sparse (data, pos) 
  1280.      DCTBLOCK data;
  1281.      int pos;
  1282. {
  1283.   j_rev_dct(data);
  1284. }
  1285.  
  1286. void
  1287. j_rev_dct (data)
  1288.   DCTBLOCK data;
  1289. {
  1290.   INT32 tmp0, tmp1, tmp2, tmp3;
  1291.   INT32 tmp10, tmp11, tmp12, tmp13;
  1292.   INT32 z1, z2, z3, z4, z5;
  1293.   register DCTELEM *dataptr;
  1294.   int rowctr;
  1295.   SHIFT_TEMPS
  1296.  
  1297.   /* Pass 1: process rows. */
  1298.   /* Note results are scaled up by sqrt(8) compared to a true IDCT; */
  1299.   /* furthermore, we scale the results by 2**PASS1_BITS. */
  1300.  
  1301.   dataptr = data;
  1302.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  1303.     /* Due to quantization, we will usually find that many of the input
  1304.      * coefficients are zero, especially the AC terms.  We can exploit this
  1305.      * by short-circuiting the IDCT calculation for any row in which all
  1306.      * the AC terms are zero.  In that case each output is equal to the
  1307.      * DC coefficient (with scale factor as needed).
  1308.      * With typical images and quantization tables, half or more of the
  1309.      * row DCT calculations can be simplified this way.
  1310.      */
  1311.  
  1312.     if ((dataptr[1] | dataptr[2] | dataptr[3] | dataptr[4] |
  1313.      dataptr[5] | dataptr[6] | dataptr[7]) == 0) {
  1314.       /* AC terms all zero */
  1315.       DCTELEM dcval = (DCTELEM) (dataptr[0] << PASS1_BITS);
  1316.       
  1317.       dataptr[0] = dcval;
  1318.       dataptr[1] = dcval;
  1319.       dataptr[2] = dcval;
  1320.       dataptr[3] = dcval;
  1321.       dataptr[4] = dcval;
  1322.       dataptr[5] = dcval;
  1323.       dataptr[6] = dcval;
  1324.       dataptr[7] = dcval;
  1325.       
  1326.       dataptr += DCTSIZE;    /* advance pointer to next row */
  1327.       continue;
  1328.     }
  1329.  
  1330.     /* Even part: reverse the even part of the forward DCT. */
  1331.     /* The rotator is sqrt(2)*c(-6). */
  1332.  
  1333.     z2 = (INT32) dataptr[2];
  1334.     z3 = (INT32) dataptr[6];
  1335.  
  1336.     z1 = MULTIPLY(z2 + z3, FIX(0.541196100));
  1337.     tmp2 = z1 + MULTIPLY(z3, - FIX(1.847759065));
  1338.     tmp3 = z1 + MULTIPLY(z2, FIX(0.765366865));
  1339.  
  1340.     tmp0 = ((INT32) dataptr[0] + (INT32) dataptr[4]) << CONST_BITS;
  1341.     tmp1 = ((INT32) dataptr[0] - (INT32) dataptr[4]) << CONST_BITS;
  1342.  
  1343.     tmp10 = tmp0 + tmp3;
  1344.     tmp13 = tmp0 - tmp3;
  1345.     tmp11 = tmp1 + tmp2;
  1346.     tmp12 = tmp1 - tmp2;
  1347.     
  1348.     /* Odd part per figure 8; the matrix is unitary and hence its
  1349.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  1350.      */
  1351.  
  1352.     tmp0 = (INT32) dataptr[7];
  1353.     tmp1 = (INT32) dataptr[5];
  1354.     tmp2 = (INT32) dataptr[3];
  1355.     tmp3 = (INT32) dataptr[1];
  1356.  
  1357.     z1 = tmp0 + tmp3;
  1358.     z2 = tmp1 + tmp2;
  1359.     z3 = tmp0 + tmp2;
  1360.     z4 = tmp1 + tmp3;
  1361.     z5 = MULTIPLY(z3 + z4, FIX(1.175875602)); /* sqrt(2) * c3 */
  1362.     
  1363.     tmp0 = MULTIPLY(tmp0, FIX(0.298631336)); /* sqrt(2) * (-c1+c3+c5-c7) */
  1364.     tmp1 = MULTIPLY(tmp1, FIX(2.053119869)); /* sqrt(2) * ( c1+c3-c5+c7) */
  1365.     tmp2 = MULTIPLY(tmp2, FIX(3.072711026)); /* sqrt(2) * ( c1+c3+c5-c7) */
  1366.     tmp3 = MULTIPLY(tmp3, FIX(1.501321110)); /* sqrt(2) * ( c1+c3-c5-c7) */
  1367.     z1 = MULTIPLY(z1, - FIX(0.899976223)); /* sqrt(2) * (c7-c3) */
  1368.     z2 = MULTIPLY(z2, - FIX(2.562915447)); /* sqrt(2) * (-c1-c3) */
  1369.     z3 = MULTIPLY(z3, - FIX(1.961570560)); /* sqrt(2) * (-c3-c5) */
  1370.     z4 = MULTIPLY(z4, - FIX(0.390180644)); /* sqrt(2) * (c5-c3) */
  1371.     
  1372.     z3 += z5;
  1373.     z4 += z5;
  1374.     
  1375.     tmp0 += z1 + z3;
  1376.     tmp1 += z2 + z4;
  1377.     tmp2 += z2 + z3;
  1378.     tmp3 += z1 + z4;
  1379.  
  1380.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  1381.  
  1382.     dataptr[0] = (DCTELEM) DESCALE(tmp10 + tmp3, CONST_BITS-PASS1_BITS);
  1383.     dataptr[7] = (DCTELEM) DESCALE(tmp10 - tmp3, CONST_BITS-PASS1_BITS);
  1384.     dataptr[1] = (DCTELEM) DESCALE(tmp11 + tmp2, CONST_BITS-PASS1_BITS);
  1385.     dataptr[6] = (DCTELEM) DESCALE(tmp11 - tmp2, CONST_BITS-PASS1_BITS);
  1386.     dataptr[2] = (DCTELEM) DESCALE(tmp12 + tmp1, CONST_BITS-PASS1_BITS);
  1387.     dataptr[5] = (DCTELEM) DESCALE(tmp12 - tmp1, CONST_BITS-PASS1_BITS);
  1388.     dataptr[3] = (DCTELEM) DESCALE(tmp13 + tmp0, CONST_BITS-PASS1_BITS);
  1389.     dataptr[4] = (DCTELEM) DESCALE(tmp13 - tmp0, CONST_BITS-PASS1_BITS);
  1390.  
  1391.     dataptr += DCTSIZE;        /* advance pointer to next row */
  1392.   }
  1393.  
  1394.   /* Pass 2: process columns. */
  1395.   /* Note that we must descale the results by a factor of 8 == 2**3, */
  1396.   /* and also undo the PASS1_BITS scaling. */
  1397.  
  1398.   dataptr = data;
  1399.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  1400.     /* Columns of zeroes can be exploited in the same way as we did with rows.
  1401.      * However, the row calculation has created many nonzero AC terms, so the
  1402.      * simplification applies less often (typically 5% to 10% of the time).
  1403.      * On machines with very fast multiplication, it's possible that the
  1404.      * test takes more time than it's worth.  In that case this section
  1405.      * may be commented out.
  1406.      */
  1407.  
  1408. #ifndef NO_ZERO_COLUMN_TEST
  1409.     if ((dataptr[DCTSIZE*1] | dataptr[DCTSIZE*2] | dataptr[DCTSIZE*3] |
  1410.      dataptr[DCTSIZE*4] | dataptr[DCTSIZE*5] | dataptr[DCTSIZE*6] |
  1411.      dataptr[DCTSIZE*7]) == 0) {
  1412.       /* AC terms all zero */
  1413.       DCTELEM dcval = (DCTELEM) DESCALE((INT32) dataptr[0], PASS1_BITS+3);
  1414.       
  1415.       dataptr[DCTSIZE*0] = dcval;
  1416.       dataptr[DCTSIZE*1] = dcval;
  1417.       dataptr[DCTSIZE*2] = dcval;
  1418.       dataptr[DCTSIZE*3] = dcval;
  1419.       dataptr[DCTSIZE*4] = dcval;
  1420.       dataptr[DCTSIZE*5] = dcval;
  1421.       dataptr[DCTSIZE*6] = dcval;
  1422.       dataptr[DCTSIZE*7] = dcval;
  1423.       
  1424.       dataptr++;        /* advance pointer to next column */
  1425.       continue;
  1426.     }
  1427. #endif
  1428.  
  1429.     /* Even part: reverse the even part of the forward DCT. */
  1430.     /* The rotator is sqrt(2)*c(-6). */
  1431.  
  1432.     z2 = (INT32) dataptr[DCTSIZE*2];
  1433.     z3 = (INT32) dataptr[DCTSIZE*6];
  1434.  
  1435.     z1 = MULTIPLY(z2 + z3, FIX(0.541196100));
  1436.     tmp2 = z1 + MULTIPLY(z3, - FIX(1.847759065));
  1437.     tmp3 = z1 + MULTIPLY(z2, FIX(0.765366865));
  1438.  
  1439.     tmp0 = ((INT32) dataptr[DCTSIZE*0] + (INT32) dataptr[DCTSIZE*4]) << CONST_BITS;
  1440.     tmp1 = ((INT32) dataptr[DCTSIZE*0] - (INT32) dataptr[DCTSIZE*4]) << CONST_BITS;
  1441.  
  1442.     tmp10 = tmp0 + tmp3;
  1443.     tmp13 = tmp0 - tmp3;
  1444.     tmp11 = tmp1 + tmp2;
  1445.     tmp12 = tmp1 - tmp2;
  1446.     
  1447.     /* Odd part per figure 8; the matrix is unitary and hence its
  1448.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  1449.      */
  1450.  
  1451.     tmp0 = (INT32) dataptr[DCTSIZE*7];
  1452.     tmp1 = (INT32) dataptr[DCTSIZE*5];
  1453.     tmp2 = (INT32) dataptr[DCTSIZE*3];
  1454.     tmp3 = (INT32) dataptr[DCTSIZE*1];
  1455.  
  1456.     z1 = tmp0 + tmp3;
  1457.     z2 = tmp1 + tmp2;
  1458.     z3 = tmp0 + tmp2;
  1459.     z4 = tmp1 + tmp3;
  1460.     z5 = MULTIPLY(z3 + z4, FIX(1.175875602)); /* sqrt(2) * c3 */
  1461.     
  1462.     tmp0 = MULTIPLY(tmp0, FIX(0.298631336)); /* sqrt(2) * (-c1+c3+c5-c7) */
  1463.     tmp1 = MULTIPLY(tmp1, FIX(2.053119869)); /* sqrt(2) * ( c1+c3-c5+c7) */
  1464.     tmp2 = MULTIPLY(tmp2, FIX(3.072711026)); /* sqrt(2) * ( c1+c3+c5-c7) */
  1465.     tmp3 = MULTIPLY(tmp3, FIX(1.501321110)); /* sqrt(2) * ( c1+c3-c5-c7) */
  1466.     z1 = MULTIPLY(z1, - FIX(0.899976223)); /* sqrt(2) * (c7-c3) */
  1467.     z2 = MULTIPLY(z2, - FIX(2.562915447)); /* sqrt(2) * (-c1-c3) */
  1468.     z3 = MULTIPLY(z3, - FIX(1.961570560)); /* sqrt(2) * (-c3-c5) */
  1469.     z4 = MULTIPLY(z4, - FIX(0.390180644)); /* sqrt(2) * (c5-c3) */
  1470.     
  1471.     z3 += z5;
  1472.     z4 += z5;
  1473.     
  1474.     tmp0 += z1 + z3;
  1475.     tmp1 += z2 + z4;
  1476.     tmp2 += z2 + z3;
  1477.     tmp3 += z1 + z4;
  1478.  
  1479.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  1480.  
  1481.     dataptr[DCTSIZE*0] = (DCTELEM) DESCALE(tmp10 + tmp3,
  1482.                        CONST_BITS+PASS1_BITS+3);
  1483.     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp10 - tmp3,
  1484.                        CONST_BITS+PASS1_BITS+3);
  1485.     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp11 + tmp2,
  1486.                        CONST_BITS+PASS1_BITS+3);
  1487.     dataptr[DCTSIZE*6] = (DCTELEM) DESCALE(tmp11 - tmp2,
  1488.                        CONST_BITS+PASS1_BITS+3);
  1489.     dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(tmp12 + tmp1,
  1490.                        CONST_BITS+PASS1_BITS+3);
  1491.     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp12 - tmp1,
  1492.                        CONST_BITS+PASS1_BITS+3);
  1493.     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp13 + tmp0,
  1494.                        CONST_BITS+PASS1_BITS+3);
  1495.     dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(tmp13 - tmp0,
  1496.                        CONST_BITS+PASS1_BITS+3);
  1497.     
  1498.     dataptr++;            /* advance pointer to next column */
  1499.   }
  1500. }
  1501.  
  1502.  
  1503. #endif
  1504.